library(BiocManager)
BiocManager::install('sigminer')

library(VariantAnnotation)
browseVignettes('VariantAnnotation')
#VariantAnnotation
fl <- system.file("extdata", "chr22.vcf.gz", package="VariantAnnotation")
vcf <- readVcf(fl, "hg19")
vcf
header(vcf)
geno(header(vcf))
ref(vcf)[1:5]
alt(vcf)[1:5]
sapply(geno(vcf), class)
geno(header(vcf))["DS",]
DS <-geno(vcf)$DS
dim(DS)
fivenum(DS)
length(which(DS==0))/length(DS)
hist(DS[DS != 0], breaks=seq(0, 2, by=0.05),
     main="DS non-zero values", xlab="DS")
info(vcf)[1:4, 1:5]
#maftools
library(maftools)

#path to TCGA LAML MAF file
laml.maf = system.file('extdata', 'tcga_laml.maf.gz', package = 'maftools') 
#clinical information containing survival information and histology. This is optional
laml.clin = system.file('extdata', 'tcga_laml_annot.tsv', package = 'maftools') 
laml = read.maf(maf = laml.maf, clinicalData = laml.clin)
#Shows sample summry.
getSampleSummary(laml)
#Shows gene summary.
getGeneSummary(laml)
#shows clinical data associated with samples
getClinicalData(laml)
#Shows all fields in MAF
getFields(laml)
#Writes maf summary to an output file with basename laml.
write.mafSummary(maf = laml, basename = 'laml')
plotmafSummary(maf = laml, rmOutlier = TRUE, addStat = 'median', dashboard = TRUE, titvRaw = FALSE)
laml.titv = titv(maf = laml, plot = FALSE, useSyn = TRUE)
#plot titv summary
plotTiTv(res = laml.titv)
laml.seg <- system.file("extdata", "LAML_CBS_segments.tsv.gz", package = "maftools")
segSummarize_results = segSummarize(seg = laml.seg)
tcga.ab.009.seg <- system.file("extdata", "TCGA.AB.3009.hg19.seg.txt", package = "maftools")
plotCBSsegments(cbsFile = tcga.ab.009.seg)

laml.sig = oncodrive(maf = laml, AACol = 'Protein_Change', minMut = 5, pvalMethod = 'zscore')
head(laml.sig)
plotOncodrive(res = laml.sig, fdrCutOff = 0.1, useFraction = TRUE, labelSize = 0.5)

laml.pfam = pfamDomains(maf = laml, AACol = 'Protein_Change', top = 10)

mafSurvival(maf = laml, genes = 'DNMT3A', time = 'days_to_last_followup', Status = 'Overall_Survival_Status', isTCGA = TRUE)
#sigminer
library(sigminer)
browseVignettes('sigminer')
brca <- readRDS("./基因突变/brca.rds")


